% PRODUCEFIG2.M
%
% Produces figure 2 in the paper.

clear all;

param;

imp = Omega;

up_to = 32;
time = 1:up_to;

%% Response under Initial Structure 
% Sets up the Stationary version of Ireland's (2007) New Keynesian model that we use for the applications in our paper. 
%
%                        1       2     3     4         5              6          7       8        9       10 
% Variables in Y(t) = [ yhat(t) pi(t) r(t) ghat(t) E_t(yhat(t+1)) E_t(pi(t+1)) ahat(t) ehat(t) ln(z(t)) zhat(t)]'
% 
%
%                           1      2    3     4
% Variables in eps(t) = [ ea(t) ee(t) ez(t) er(t) ]'
%
%
% Variables in eta(t) = [etay(t) etapi(t)]'

% Load Dimensions
n1 = 6;
n2 = 2;
nleads = 1;
s = nleads;
nfcast = 2;
nshocks = 4;
nvars = n1+n2+nfcast;
n = nvars;
k = nfcast;

% Load Initial Matrices
% Note: row indicates equation number and column indicates variables number
GAM0 = zeros(nvars);
GAM0(1,1) = 1; 
GAM0(1,3) = 1/sigma;
GAM0(1,5) = -1;
GAM0(1,6) = -1/sigma; 
GAM0(1,7) = -(1-rho_a)/sigma;
GAM0(2,1) = -psi*sigma;
GAM0(2,2) = 1+beta*alpha;
GAM0(2,6) = -beta;
GAM0(2,8) = 1;
GAM0(2,10) = psi;
GAM0(3,1) = -rho_y;
GAM0(3,2) = -rho_pi;
GAM0(3,3) = 1;
GAM0(3,4) = -rho_g;
GAM0(4,7) = 1;
GAM0(5,8) = 1;
GAM0(6,9) = 1;
GAM0(7,1) = -1;
GAM0(7,4) = 1;
GAM0(8,9) = -1;
GAM0(8,10) = 1;
GAM0(9,1) = 1;
GAM0(10,2) = 1;

GAM1 = zeros(nvars);
GAM1(2,2) = alpha;
GAM1(3,3) = rho_r;
GAM1(4,7) = rho_a;
GAM1(5,8) = rho_e;
GAM1(6,9) = rho_z;
GAM1(7,1) = -1;
GAM1(9,5) = 1;
GAM1(10,6) = 1;

C = zeros(nvars,1);
C(1,1) = -(1/sigma)*log(beta);
C(2,1) = (1+beta*alpha-alpha-beta)*pitarget;
C(3,1) = (1-rho_r)*(pitarget-log(beta)) - rho_pi*pitarget;
C(6,1) = (1-rho_z)*log(z);
C(8,1) = -log(z);

PSI = zeros(nvars,nshocks);
PSI(3,4) = 1;
PSI(4,1) = 1;
PSI(5,2) = 1;
PSI(6,3) = 1;

PPI = zeros(nvars,nfcast);
PPI(9,1) = 1;
PPI(10,2) = 1;

[S0_i, S1_i, S2_i, S3_i, unique_i] = Smats(C, GAM0, GAM1, PSI, PPI);
y0 = inv(eye(n) - S1_i)*S0_i;

yirfini = zeros(n,up_to);

% Impulse response to ea(t)
yirfini(:,1) = inv(eye(n)-S1_i)*S0_i + S2_i*imp(:,1);
for t = 1:(up_to-1)
    yirfini(:,t+1) = S0_i + S1_i*yirfini(:,t);
end

%% Response under Final Structure 
GAM0n = zeros(nvars);
GAM0n(1,1) = 1; 
GAM0n(1,3) = 1/sigma_n;
GAM0n(1,5) = -1;
GAM0n(1,6) = -1/sigma_n; 
GAM0n(1,7) = -(1-rho_a_n)/sigma_n;
GAM0n(2,1) = -psi_n*sigma_n;
GAM0n(2,2) = 1+beta_n*alpha_n;
GAM0n(2,6) = -beta_n;
GAM0n(2,8) = 1;
GAM0n(2,10) = psi_n;
GAM0n(3,1) = -rho_y_n;
GAM0n(3,2) = -rho_pi_n;
GAM0n(3,3) = 1;
GAM0n(3,4) = -rho_g_n;
GAM0n(4,7) = 1;
GAM0n(5,8) = 1;
GAM0n(6,9) = 1;
GAM0n(7,1) = -1;
GAM0n(7,4) = 1;
GAM0n(8,9) = -1;
GAM0n(8,10) = 1;
GAM0n(9,1) = 1;
GAM0n(10,2) = 1;

GAM1n = zeros(nvars);
GAM1n(2,2) = alpha_n;
GAM1n(3,3) = rho_r_n;
GAM1n(4,7) = rho_a_n;
GAM1n(5,8) = rho_e_n;
GAM1n(6,9) = rho_z_n;
GAM1n(7,1) = -1;
GAM1n(9,5) = 1;
GAM1n(10,6) = 1;

Cn = zeros(nvars,1);
Cn(1,1) = -(1/sigma_n)*log(beta_n);
Cn(2,1) = (1+beta_n*alpha_n-alpha_n-beta_n)*pitarget_n;
Cn(3,1) = (1-rho_r_n)*(pitarget_n-log(beta_n)) - rho_pi_n*pitarget_n;
Cn(6,1) = (1-rho_z_n)*log(z_n);
Cn(8,1) = -log(z_n);

PSIn = zeros(nvars,nshocks);
PSIn(3,4) = 1;
PSIn(4,1) = 1;
PSIn(5,2) = 1;
PSIn(6,3) = 1;

PPIn = zeros(nvars,nfcast);
PPIn(9,1) = 1;
PPIn(10,2) = 1;

[n l] = size(PSI);

% Reduced form solution matrices of final structure.
[S0_f, S1_f, S2_f, S3_f, unique_f] = Smats(Cn, GAM0n, GAM1n, PSIn, PPIn);
y0_new = inv(eye(n) - S1_f)*S0_f;

% Impulse response to ea(t)
yirffin(:,1) = inv(eye(n)-S1_f)*S0_f + S2_f*imp(:,1);
for t = 1:(up_to-1)
    yirffin(:,t+1) = S0_f + S1_f*yirffin(:,t);
end

%% Response with T = 4
% Compute transition period
% Store endogenous variables in the transition period
y_transition_4 = zeros(n,T);
% Load unanticipated shocks (use in stochastic simulations)
et = zeros(l,2*T); % size
% Store variables for overall response
yirfT4 = zeros(n,up_to);

% Construct sequences of time-varying matrices 
GAM0s = zeros(n*(T+1),n);
GAM1s = zeros(n*(T+1),n);
Cs = zeros(n*(T+1),1);
PSIs = zeros(n*(T+1),l);
PPIs = zeros(n*(T+1),k);

GAM0s(1:n*T,:) = repmat(GAM0,T,1);
GAM0s(n*T+1:end,:) = GAM0n;

GAM1s(1:n*T,:) = repmat(GAM1,T,1);
GAM1s(n*T+1:end,:) = GAM1n;

Cs(1:n*T,:) = repmat(C,T,1);
Cs(n*T+1:end,:) = Cn;

PSIs(1:n*T,:) = repmat(PSI,T,1);
PSIs(n*T+1:end,:) = PSIn;

PPIs(1:n*T,:) = repmat(PPI,T,1);
PPIs(n*T+1:end,:) = PPI;

complex_tol = 1e-6;
eps_a = zeros(l,T); % Load anticipated shocks

% IRF up to period ta-1
yirfT4(:,1:ta-1) = yirfini(:,1:ta-1);

% For initial period of the transition
[yseq, S0, S1, S2, Sy, Sf, Sz, A, B] = LREAnt(n1,n2,s,T,GAM0s,GAM1s,Cs,PSIs,PPIs,yirfT4(:,ta-1),eps_a, et(:,ta),complex_tol);
y0 = yseq(:,1);
y_transition_4(:,1) = yseq(:,1);
% Loop is useful for stochastic simulation.
for t = 2:T
    eps_a = zeros(l,T-t+1);
    GAM0s(1:n,:) = [];
    GAM1s(1:n,:) = [];
    Cs(1:n,:) = [];
    PSIs(1:n,:) = [];
    PPIs(1:n,:) = [];
   [yseq, S0, S1, S2, Sy, Sf, Sz, A, B] = LREAnt(n1,n2,s,T-t+1,GAM0s,GAM1s,Cs,PSIs,PPIs,y0,eps_a, et(:,ta+t-1),complex_tol);
    y0 = yseq(:,1);
    y_transition_4(:,t) = yseq(:,1);
end

yirfT4(:,ta:T+ta-1) = y_transition_4;

% End of Response with final structure

for t = T+ta:up_to
     yirfT4(:,t) = S0_f + S1_f * yirfT4(:,t-1);
end


%% Response with T = 24
T = 20;
% Compute transition period
% Store endogenous variables in the transition period
y_transition_24 = zeros(n,T);
% Load unanticipated shocks (use in stochastic simulations)
et = zeros(l,2*T); % size
% Store variables for overall response
yirfT24 = zeros(n,up_to);

% Construct sequences of time-varying matrices 
GAM0s = zeros(n*(T+1),n);
GAM1s = zeros(n*(T+1),n);
Cs = zeros(n*(T+1),1);
PSIs = zeros(n*(T+1),l);
PPIs = zeros(n*(T+1),k);

GAM0s(1:n*T,:) = repmat(GAM0,T,1);
GAM0s(n*T+1:end,:) = GAM0n;

GAM1s(1:n*T,:) = repmat(GAM1,T,1);
GAM1s(n*T+1:end,:) = GAM1n;

Cs(1:n*T,:) = repmat(C,T,1);
Cs(n*T+1:end,:) = Cn;

PSIs(1:n*T,:) = repmat(PSI,T,1);
PSIs(n*T+1:end,:) = PSIn;

PPIs(1:n*T,:) = repmat(PPI,T,1);
PPIs(n*T+1:end,:) = PPI;

complex_tol = 1e-6;
eps_a = zeros(l,T); % Load anticipated shocks

% IRF up to period ta-1
yirfT24(:,1:ta-1) = yirfini(:,1:ta-1);

% For initial period of the transition
[yseq, S0, S1, S2, Sy, Sf, Sz, A, B] = LREAnt(n1,n2,s,T,GAM0s,GAM1s,Cs,PSIs,PPIs,yirfT24(:,ta-1),eps_a, et(:,ta),complex_tol);
y0 = yseq(:,1);
y_transition_24(:,1) = yseq(:,1);
% Loop is useful for stochastic simulation.
for t = 2:T
    eps_a = zeros(l,T-t+1);
    GAM0s(1:n,:) = [];
    GAM1s(1:n,:) = [];
    Cs(1:n,:) = [];
    PSIs(1:n,:) = [];
    PPIs(1:n,:) = [];
   [yseq, S0, S1, S2, Sy, Sf, Sz, A, B] = LREAnt(n1,n2,s,T-t+1,GAM0s,GAM1s,Cs,PSIs,PPIs,y0,eps_a, et(:,ta+t-1),complex_tol);
    y0 = yseq(:,1);
    y_transition_24(:,t) = yseq(:,1);
end

yirfT24(:,ta:T+ta-1) = y_transition_24;

% End of Response with final structure

for t = T+ta:up_to
     yirfT24(:,t) = S0_f + S1_f * yirfT24(:,t-1);
end


%% Plot Figure 2
% Transform Units
inflation_ini = yirfini(2,:)*400;
inflation_fin = yirffin(2,:)*400;
inflation_T4 = yirfT4(2,:)*400;
inflation_T24 = yirfT24(2,:)*400;
figure('Name','Figure 2 in the paper: Demand shock','NumberTitle','off')
plot(time,inflation_ini,'b--',time,inflation_fin,'r--',time,inflation_T4,'g',time,inflation_T24,'k--','LineWidth',lw),
         title('Inflation','fontsize',14), ylabel('%'),legend('Initial','Final','Path4','Path24'),grid,set(gca,'xtick',[ta, 4+ta, 20+ta])





